Week 6: Visualisation of Spatial Data – Maps

1 Introduction

• In this Lab, we will learn to draw some maps.

2 USA map

R package maps and mapdata provide some maps. We will use them for the USA map

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.6.2

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(maps)
library(mapdata)
library(sp)
usa <- map_data("usa")
dim(usa)

## [1] 7243    6

head(usa)

##        long      lat group order region subregion
## 1 -101.4078 29.74224     1     1   main      <NA>
## 2 -101.3906 29.74224     1     2   main      <NA>
## 3 -101.3620 29.65056     1     3   main      <NA>
## 4 -101.3505 29.63911     1     4   main      <NA>
## 5 -101.3219 29.63338     1     5   main      <NA>
## 6 -101.3047 29.64484     1     6   main      <NA>

tail(usa)

##           long      lat group order         region subregion
## 7247 -122.6187 48.37482    10  7247 whidbey island      <NA>
## 7248 -122.6359 48.35764    10  7248 whidbey island      <NA>
## 7249 -122.6703 48.31180    10  7249 whidbey island      <NA>
## 7250 -122.7218 48.23732    10  7250 whidbey island      <NA>
## 7251 -122.7104 48.21440    10  7251 whidbey island      <NA>
## 7252 -122.6703 48.17429    10  7252 whidbey island      <NA>

First, we plot usa map

library(gridExtra)
gg0 = ggplot() + geom_polygon(data = usa, aes(x=long, y = lat, group = group)) + 
  coord_fixed(1.3)
gg1 <- ggplot() + 
  geom_polygon(data = usa, aes(x=long, y = lat, group = group), fill = "violet", color = "blue") + 
  coord_fixed(1.3)
labs <- data.frame(
  long = c(-122.064873, -122.306417),
  lat = c(36.951968, 47.644855),
  names = c("SWFSC-FED", "NWFSC"),
  stringsAsFactors = FALSE
)  
gg2 = gg1 + 
  geom_point(data = labs, aes(x = long, y = lat), color = "yellow", size = 2)

grid.arrange(gg0, gg1, gg2, nrow = 1)

• Key points
  • “group = group” is essential, you can try code without the group argument
  • coord_fixed(): It fixes the relationship between one unit in the y direction and one unit in the x direction.

2.1 US state

### Draw plots for all USA States
states <- map_data("state")
ggplot(data = states) + 
  geom_polygon(aes(x = long, y = lat, fill = region, group = group), color = "white") + 
  coord_fixed(1.3) + guides(fill=FALSE) 

### Draw a subset of USA states
### Use function subset to select states. 
west_coast <- subset(states, region %in% c("california", "oregon", "washington"))
ggplot(data = west_coast) + 
  geom_polygon(aes(x = long, y = lat, group = group), fill = "palegreen", color = "black") + 
  coord_fixed(1.3)

Can you tell the difference between color and fill argument?

3 Victoria maps

• Most shapefile can found in https://data.gov.au.
• All LGA: VICLGA folder – 92 elements
• All suburb: VICSub folder – 2973 elements
• LGA in Greater Melbourne – 32 element (1st element is greater melbourne)

All datasets are read into as SpatialPolygonsDataFrame

library(rgdal)
library(sp)
VICLGA = readOGR(dsn = "datasets/VICLGA/VIC_LOCALITY_POLYGON_shp.shp", "VIC_LOCALITY_POLYGON_shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/tingjinc/Dropbox/2_MAST90122/90122LabGit/datasets/VICLGA/VIC_LOCALITY_POLYGON_shp.shp", layer: "VIC_LOCALITY_POLYGON_shp"
## with 92 features
## It has 10 fields

VICSub = readOGR(dsn="datasets/VICSub/VIC_LOCALITY_POLYGON_shp.shp", "VIC_LOCALITY_POLYGON_shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/tingjinc/Dropbox/2_MAST90122/90122LabGit/datasets/VICSub/VIC_LOCALITY_POLYGON_shp.shp", layer: "VIC_LOCALITY_POLYGON_shp"
## with 2973 features
## It has 12 fields

MELLGA = readOGR(dsn = "datasets/Mel/Aus-Melbourne02.shp", "Aus-Melbourne02")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/tingjinc/Dropbox/2_MAST90122/90122LabGit/datasets/Mel/Aus-Melbourne02.shp", layer: "Aus-Melbourne02"
## with 32 features
## It has 14 fields

• You will cut SpatialPolygons you wants – Essentially manipulation of SpatialPolygons.
• One useful funciton in crop in package raster
  • crop(x, y)
• crop returns a geographic subset of an object x as specified by an Extent object y
  • ?crop for more details
• Here we will create
• Task 1: all surburbs in Greater Melbourne
• Task 2: all surburbs in City of Melbourne

library(raster)
library(mapview)

## Warning: package 'mapview' was built under R version 3.6.2

# Before crop, first check CRS
VICLGA@proj4string

## CRS arguments:
##  +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

VICSub@proj4string

## CRS arguments: +proj=longlat +ellps=GRS80 +no_defs

MELLGA@proj4string

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

# Making all CRS the same
crs0 = VICLGA@proj4string
MELLGA = spTransform(MELLGA, crs0)
MELLGA@proj4string

## CRS arguments:
##  +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

## Task 1
GMEL = MELLGA[1, ]   # Greater Melbourne
GMELSub = crop(VICSub, GMEL)
mapView(GMELSub)      

## Task 2: this code takes a while to run. 
i = which(MELLGA$name == "City of Melbourne")
CMEL = MELLGA[i, ]
#CMELSub = crop(VICSub, CMEL)      
#mapView(CMELSub) 

• There are more shapefiles in PlanMelbourne Folder, which can be downloaded from https://www.planmelbourne.vic.gov.au/maps/spatial-data
• I have not explore all of shapefiles
• It is quite useful.

• The following gives regions of Victoria and CBD

  #SReg = readOGR(dsn="datasets/PlanMelbourne/Administrative/State Regions_region.shp", "State Regions_region")
  
  #library(mapview)
  #mapView(SReg)
  #CBD = readOGR(dsn="2_ggplot/PlanMelbourne/Central City Urban Renewal and Precincts/Central Business District_region.shp", "Central Business District_region")
  #mapView(CBD)